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ABSTRACT 


The hydrodynamic problem for multi-body systems in waves has complexity in both wave structure 
interaction and coupled dynamics with interconnection. A consistent approach is developed for calculating the 
dynamic response of multi-body floating systems. This method is based on multi-rigid-body hydrodynamics and 
Euler-beam bending theory and can be used to deal with rigid or flexible multi-body floating systems in both 
frequency and time domain. Within the framework of this consistent approach, a stiffness approximation method 
is proposed to calculate the natural frequency and corresponding oscillation mode of a hinged multi-module 
structure. This is significant for both traditional multi-body floating system (for the purpose of structural safety) 
and wave energy converters (for the purpose of power efficiency). The hydroelasticity theory proposed by Lu et al. 
(2016), which is used to deal with continuous flexible structures with cross section of simple shapes, is extended 
to be applicable for multi-module flexible structures with cross section of arbitrary shapes and varied geometric 
features in longitudinal direction. Finally, an effective strategy is established to analyse in time domain the 
dynamic response of a hinged multi-module flexible structure with moving point loads on the upper surface in 


waves, representative of vehicles moving on a floating bridge or airport. 
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1. Introduction 


Multi-body floating system has been widely designed and used for a variety of purpose in both coastal and 
ocean region, such as permanent floating structures including multi-module floating airport and bridge, floating 
LNG terminal clusters, hinged raft-type wave energy converters, wave energy converter arrays; and temporary 
multi-body floating system including offloading operation between FPSO (or FLNG) and LNG carrier, assembly 
process of huge floating structures and a crane barge in the vicinity of an offshore platform. The hydrodynamic 
problem for multi-body systems in waves has long been studied due to its complexity lying in not only wave and 
(rigid or flexible) structure interaction but also the coupled dynamics with interconnection considered. Various 
approaches have been proposed to deal with hydrodynamic problems of (rigid or flexible) multi-body floating 


structures, in both frequency and time domain. 
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Many researchers have adopted the linear frequency domain potential flow model to deal with multi-body 
hydrodynamic problems. Earlier work on two adjacent bodies in waves without connection in-between done by 
Ohkushu (1974), Kodan (1984) and Fang and Kim (1986) highlighted the effects of multi-body hydrodynamic 
interactions. More recently, Inoue et al. (2000) and Kashiwagi et al. (2005) focused on the drift force and moment 
for the two adjacent floating bodies. For wave energy converter arrays, frequency domain method is used to 
establish the motion equations and then investigate effects of multi-body hydrodynamic interaction on the power 
capture efficiency. Relevant work was undertaken by Babarit (2010), Babarit (2013) and Borgarino et al. (2012). 
The difference between traditional multi-body hydrodynamics and multi-body hydrodynamics for wave energy 
converter arrays is that the latter needs an equivalent power-take-off damping coefficient to be added in the 
motion equation. The above-mentioned work focused mainly on freely floating multi-body structures without 
constraints in-between, which is actually the basic multi-body hydrodynamic problem. 

For interconnected multi-module floating structures, the hydrodynamic problems become more complex due 
to the coupling effects caused by the interconnection (the interconnection form may be hinges or moorings). 
Newman (1994) proposed a mode expansion method to deal with the hydrodynamic response of (rigidly) hinged 
floating bodies within the framework of linear frequency-domain method. In his study, the multi-module structure 
is considered as a whole with number of modes equal to the number of degrees of freedom (DOFs). Then the 
dynamic response of a hinged structure can be obtained based on the following two steps: (1) Evaluation of the 
hydrodynamic coefficients (added mass, radiation damping and wave excitation force) for each mode; and (2) 
Solving the coupling modal equation to obtain the hydrodynamic response of a hinged structure. This method was 
also adopted by Taghipour and Moan (2008) to investigate the performance of a multi-body wave energy 
converter in absorbing wave energy and its dynamic behaviour in multi-directional waves. Nevertheless, for 
floating systems with complex constraints or geometrically different modules, the structure cannot be regarded as 
a single body and thus the mode expansion method may be invalid. Gou et al. (2004) utilized a different approach 
by introducing a constraint matrix to study the interaction between waves and (rigidly) hinged two-module 
structure. The hydrodynamic coefficients were calculated for two modules without considering the connection, 
after which the displacement continuity condition was considered. Then the motion equation of the interconnected 
two-module structure in waves was established and the dynamic response was obtained. Their results compared 
well with those calculated by Newman (1994) except that at high frequency region, some differences were found 
between the two approaches. Gou et al. attributed this difference to the wave frequency close to the natural 
frequency of the system. Based on the similar concept, Zheng et al. (2016a, b) investigated the maximum 
theoretical power absorption of connected floating bodies without and with motion constraints. Feng and Bai 
(2017) proposed a constraint matrix to model interconnections between floating bodies. 

Although in most cases, the floating system are assumed to be rigid, structural deformation may have 
significant influences on fluid field for some floating structures such as very large floating structures (VLFSs). 
The mode expansion method, which is a similar concept shown in Newman (1994), has been adopted by many 
researchers to investigate the hydroelastic response of multi-module flexible structures. In Newman’s work, the 
structure is assumed to be rigid and only rigid motion mode was considered. However, for flexible structures, all 
the modes induced by both rigid body motion and structural deformation are taken into account. Relevant 
researches have been carried out by Riggs et al. (2000), Malenica et al. (2003), Fu et al. (2007), Gao et al. (2011) 
and Sengupta et al. (2017). Unlike mode expansion scheme, Lu et al. (2016) proposed a new method to deal with 


hydroelastic response of a continuous flexible structure by combining multi-rigid-body hydrodynamics and 
Euler-Bemoulli beam assumption. This approach was adopted by Xu et al. (2017) together with constraint 
condition formulated by Gou et al. (2004) to study the hydroelastic response of (rigidly) hinged VLFSs. Good 
agreement was obtained with results given by Fu et al. (2007). 

Frequency domain method is convenient and remains dominantly utilized in offshore industry. But it is 
limited to the framework of linear assumption. If nonlinear external loading introduced by mooring, fender, riser 
or hydraulic power-take-off system is considered, time domain approach is an appropriate solution. Although fully 
nonlinear potential flow model or viscous numerical wave tank (NWT) technique is capable of calculating wave 
fields and body motions simultaneously in time domain, hybrid frequency/time domain method (i.e. the time 
domain method based on impulse response function (IRF)), which was first introduced by Cummins (1962), is 
still favoured due to its convenience and low computational cost. Using IRF based time domain approach, Koo 
and Kim (2005) simulated the side-by-side offloading operation of two platforms with mooring, riser and fender 
considered. Results show that the coupling effects of two vessels significantly influences the sway and roll 
motions. Hong et al. (2005) adopted a generalized mode approach to analyse the motion and drift force of 
side-by-side moored multiple vessels in time domain. The concept of generalized mode is extended by Tuitman et 
al. (2012) for time-domain seakeeping computations including the dynamic response of flexible barges. O’Cathain 
et al. (2008) utilized Newton-Euler equations with eliminated constraints to capture the dynamics of a 
hinged-barge wave energy device in time domain. Zheng et al. (2015) established the motion equation of two-raft 
wave energy converter in time domain using hybrid frequency/time domain approach together with a constraint 
matrix. Then the effects of Coulomb damping, radius of gyration and latching control on power efficiency were 
discussed in detail. 

The aim of this study is to develop a consistent approach capable of calculating the dynamic response of 
multi-body floating systems based on previous work by Lu et al. (2016) and Xu et al. (2017). This method can be 
used to deal with rigid or flexible multi-body floating systems (each module is freely floating or interconnected) 
in both frequency and time domain. In the present work, three main contributions are made, including (a) 
proposing an approach (stiffness approximation method, SAM) for calculating the natural frequency and the 
corresponding motion mode of interconnected multi-body floating system. The determination of natural frequency 
is crucial to not only traditional multi-body system (to reduce response by setting natural frequency outside the 
range of wave energy) but also wave energy converters (to improve power capture efficiency by tuning natural 
frequency to wave frequency); (b) extending the hydroelastic theory proposed by Lu et al. (2016) to be applicable 
for flexible structure with arbitrary shapes of cross section; and (c) investigating the hydroelastic response of 
interconnected VLFS with moving point loads on body surface in time domain, representative of vehicles moving 


on floating bridge or airport. 


2. Hydrodynamic analysis of rigid multi-body floating structures. 


2.1. Equations of motion for rigid multi-body floating structures 
The fluid around floating structures is assumed to ideal (i.e., uniform, continuous, inviscid, incompressible 
and irrotational). Then the fluid behaviour can be described by the velocity potential, øx, y, z, £) (the xoy plane is 


located on the still water surface and z is positive upwards and ¢ is time). Within the framework of linear 


assumption (i.e. the wave amplitude is small relative to wave length and body dimension), the velocity potential 


can be decomposed into three parts as follows: 


p(x, y, z, t)=6(% y, zt) +b (% Y, zt) + Ge (% ys z t) =Re{[a, + Dp +9, Je} (1) 


where, œ is the (circular) wave frequency, Ø, @ and øg denotes, respectively, the incident wave potential, 
diffraction wave potential, and radiation wave potential. ør, gp and øp are the complex amplitude of corresponding 


wave velocity potential, which satisfy the following boundary conditions: 


Vo =0, in Q 
2 0p 
-o p+ g—=0, on Sp 
OZ 
oe =0, on Sp 
Oz 
elp, +P) _ 9 és Ss, (2) 
Oz : 
OF i on S, (k =1, 2, ..., m) 
ôn, t 
Oa o on Ds S, 
on; jal jek ` 
Tee ie oa ae on S, 
r> or g 
S S 
F Module 1 Module 2 oes —-@®@-]| Module m F 
Sı S2 Sin | 


connection 1 connection 2 connection m-1 


ae 


Fig. 1. Definition of fluid and structure boundaries 


where (see Fig. 1.) Q is the fluid domain, and Sp, Sp, and S» are the free surface, bottom surface, and the boundary 
surface at infinity of the fluid, respectively. Sx (Sj) represents the wetted body surface of the k™ g™) module of 


multi-body floating system (k, j=1, 2, ... , m; jk). fig represents the outward-directed unit vector normal to the 


wetted surface of the K” module. Vs, is the velocity of a given point on the wetted surface of the K” module, and o 
is the velocity potential (In Eq.(2), g can be replaced by @, @p or Øp). r is the distance between the far-field point 
and the source point. After the velocity potential g is obtained, the added mass and the radiation damping of these 
modules as well as the wave excitation force can be calculated. 

In the frequency domain, the excitation forces are related to the incident and diffracted wave potentials as 


follows: 


F. = pial | (p, +0) ñas (3) 


So k 


where i, Foii pand So, are the imaginary unit, wave excitation forces, fluid density and average wetted surfaces, 
respectively. 

The added mass and radiation damping is given by 

i e 
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The equation of motion of freely floating multi-body system (In Fig.1, the multi-body system comprises m 

modules and each module has a six degree-of-freedom motion) in the frequency domain for a unitary wave 


amplitude and a wave frequency œ is given as follows: 
(-0 ([M, ]+ [a4 ]) -i@[bu ]+[C. ]) {u} + > > (o [ay |—ia[ by u}= {fa} E =L 2.5m) (5) 


where [M;] is the mass (or inertia moment) matrix of the k” module, [ax]is the added mass (or moment) matrix of 
the K” module caused by the motion of the module itself, [bx] is the radiation damping matrix of the K” module 
caused by the motion of the module itself, [C;] is the hydrostatic stiffness matrix of the k module, {fa} is the 
wave excitation force (or moment) of the k” module, and {u,} is the six DOF displacement of the K” module 
expressed as (Uzy Uky Uko Oy Po ny. The dimension is 6X6 for [Mi], [ar], [Duc], [Ck] [ax;] and [b;;] and 6x1 for 
{ux} and {fax} (k, j= 1, 2, ..., m). It should be noted that in Eq.(5), the constraint of displacement due to the 
existence of connections is not considered. 

If the displacement constraint due to the existence of connections is taken into account, Eq. (5) can be 


modified as follows, 
(-0 ([M, ]+ [au ]) -io[bu ]+[C. ]) {un }+ 2 È (o | ay |-ial b, » |) {u effa {Fa}. k=12,...m) (6 


where {Fi} (the dimension is 6X1) is the force and moment acting on the center of gravity of floating structure 


caused by the connection piece. 


2.2. An approach for determining natural frequency of interconnected multi -body structures 
For simplicity and clarity, in the following analysis, the number of modules for multi-body floating structures 
is restricted to 2. The results obtained are easily and straightforward to be extended to multi-body structures of 


more than 2 modules. Then for a two-module floating structure, Eq. (6) can be simplified as 
(o ete [an] J a he Pas hea fa (7) 
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As shown in Fig. 2, oxyz is earth fixed coordinate system with xoy located on the still water surface and z 


or 


being positive upwards. 01x1Z; and Ox2yxz2 are body-fixed coordinate system for module 1 and module 2, 


respectively, with o; and 02 being the center of gravity. 
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Fig. 2. A schematic of interconnected two-module floating structure 


We suppose the connection form is the common rigid hinge connection and then give a detailed explanation 
of the force {F,,} and {F,2}. The idea shown in this section can be easily extended to other forms of connection. 


At the connection, only the rotation around y axis is unconstrained. Then according to the displacement continuity 


at the connection, we have 
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Eq. (9) can be re-written in the matrix form with the introduction of constraint matrix [L], 
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If we denote the force (and moment) on the connection as {f-} (the dimension of {fc} is 5x1 as the 
moment due to the rotation around y axis is zero and removed from the force vector), then the force on the 


connection can be transferred to the center of gravity of each connected body as 


teal =[L] {fc} (11) 


The proof of Eq. (11) is given in Appendix A. Then by combining Eq. (8), (10) and (11), we have the motion 
equation for (rigidly) hinged two-module floating system, 


[K Jou. —[L} 3 es tert 
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Next we will introduce an approach to calculate the natural frequency and corresponding motion mode for 
the (rigidly) hinged two-module structure. Here we denote the method as ‘stiffness approximation method (SAM)’. 
For the hinge connection, we assume that except the free rotation around y axis, elastic deformation exists for 


other degrees of freedom (ucy Ucy Uc» Ac Yc), i.e. 
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where, Ka (s =x, y,z,@,y) 18 the stiffness coefficient; f(s =x, y, z, œ,7)is the force on the connection. The 


superscript ' and * represents, respectively, the displacement at the connection caused by the first and second 


module. Eq. (13) can be re-written as 
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In Eq. (14), the relations between displacement at the connection and at the center of gravity of the module is 
utilized, i.e. 
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If the diagonal elements of the stiffness matrix of connection approach to infinity, then the two-module floating 


structure is interconnected by a rigid hinge. Combining Eq. (11) and Eq. (14), we have the following equation, 


(ra oraa, A 


Substituting Eq. (17) into Eq. (8), we have 


(rr fi =t fo (18) 


Eq. (18) is equivalent to Eq. (12) with the diagonal elements of stiffness matrix of the connection approximating 
to infinity. Eq. (12) cannot be used to calculate the natural frequency of the (rigidly) hinged two-module structure 
due to the existence of the displacement constraint and unknown connection force. However, the natural 
frequency can be calculated if the radiation damping and wave excitation force is ignored in Eq. (18), i.e. 
[m] (EHET [Kc ILE] {a A =o u) (19) 
Mathematically speaking, the square of natural frequency is the eigenvalue of the matrix [A] and the 
corresponding oscillation mode is the eigenvector of the matrix [A]. Using Eq. (19) with diagonal elements of 
[Kc] being approximate to infinity, the natural frequency and corresponding oscillation mode can be obtained. 
Next we will use an example to validate the approach in detail. The hinged two-module floating structure in 
Newman (1994) is adopted. The structure is composed of two identical rectangular module, each of which is 40m 


long by 10m beam and 5m draft, connected end-to-end by a simple hinge with a gap of 10m between two modules 
(see Fig. 3). 
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Fig. 3. Configuration of the hinged two-module structure 


The hinge axis is located midway between two modules in the plane of the free surface. Each module is 


assumed to be in static equilibrium, with uniform mass distribution. Finally, the water depth is assumed to be 


infinite. 
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Fig. 4. Response of the hinge obtained from different approaches. (a) heave amplitude of the connector (uc,), normalized by incident 


wave amplitude (Aj); (b) Relative angle at the connection (£c), normalized by kA, (k is the wave number). 


As shown in Fig. 4, the results obtained from Eq. (12) (constraint matrix method) and stiffness approximation 
method (SAM, the stiffness value is 10'° N/m (for linear stiffness) or Nm (for angular stiffness)) are almost 
identical. It should be noted that for this case, the variation of stiffness from 10° to 10” will give almost identical 
results. Some differences exist at small wave period when the results obtained from constraint matrix method or 
SAM are compared with those obtained from mode expansion method in Newman (1994). This is due to the fact 
the natural period of the system is in the range of 5~7s and the resonance occurs at small wave period. 

We use SAM method to calculate the natural frequency and the corresponding oscillation mode of the system. 
First, we choose the wave period equal to 5s (the corresponding frequency is 1.257rad/s). And the stiffness value 
for SAM is 10” N/m (for linear stiffness) or Nm (for angular stiffness). The natural frequency obtained by solving 
Eq. (19) is œ = (0 0 O 0.6331 1.021 1.041 1.045 203.502 691.563 2929.646 5705.896 
5789.412)" (rad/s). Here the first three zero natural frequency w1,~W3 corresponds to the overall surge, sway and 
yaw motion of the system, of which the restoring force or moment is zero. The 4™ natural frequency w, is a 
complex number because the transverse metacentric height for roll motion is negative for the model. The large 
value of the 8™ to 12™ natural frequency is caused by the large stiffness value set for the corresponding degree of 
freedom (in order to approximate the real model, i.e. rigidly hinged two-module structure). Actually the natural 
frequency of interest is the 5™ to 7" one. The corresponding eigenvector (i.e. oscillation mode shape) is both listed 
in Table 1 and shown in Fig.5. It can be seen that the modal motion is the rotation around the hinge connection at 
natural frequency wW,5=1.021 rad/s, overall heave motion at w,=1.041 rad/s and overall pitch motion at 
@7=1.045 rad/s. It should be noted that the resonance does not occur as the wave frequency (1.257rad/s) is not 


equal to the natural frequency of the system. 


Table 1 Natural frequency and corresponding mode shape 


Natural frequency (rad/s) The corresponding mode shape (normalized eigenvector) 
(jx Uy Uz @ Êi Yı Ux Uy Ux, @ fo » _)' (mor rad) 
1.021 (0 O 0411 0 0.576 0 0 O 0.411 0 -0.576 0) 
1.041 © 0 0.707 0 -0.004 0 0 0 0.707 0 0.004 0)" 


1.045 (0.010 0 0.706 0 0.028 0 0.010 0 -0.706 0 0.028 0)" 
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Fig. 5. The motion mode of different natural frequencies for a two-module structure 
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Fig. 6. Determination of resonance frequency for the two-module structure 


In order to obtain wave frequency corresponding to the resonance condition of the system, we set a series of 
input wave frequency and then calculate the corresponding three natural frequencies of the floating system. The 
results are given in Fig. 6. The resonance frequency for three modal motion (i.e. rotation around the hinge, overall 
heave and pitch) is 1.014, 1.055 and 1.044 rad/s, respectively. The corresponding resonance period is within 


5.9~6.2 s. This confirms that the oscillation at low wave period shown in Fig.4 is relevant to the occurrence of 


resonance condition, i.e. wave frequency close to the natural frequency of the system. Besides, the 
above-mentioned calculation procedure can be applied to determination of the resonance condition of a hinged 
raft-type wave energy converter, which is significant for improving the power capture efficiency and considering 


the structural integrity of the system. 


3. Hydrodynamic analysis of flexible multi-body floating structures 


3.1. Equations of motion 

The hydroelasticity theory proposed by Lu et al. (2016) is adopted to investigate the hydroelastic behaviour 
of interconnected multi-body flexible structures. Compared with traditional mode expansion method (Wu, 1984), 
Lu’s approach, which is based on multi-rigid-body dynamics and beam bending, is easier to implement by just 
adding a structural stiffness matrix (taking into account the elasticity of the structure) in motion equations of rigid 
multi-body structures (for example, Eq. (6)). In this sense, Lu’s approach is a consistent theory, capable of dealing 
with dynamic response of multi-body floating system from rigid to flexible structures by turning off/on the 
structural stiffness matrix in multi-rigid-body motion equations. 

For simplicity, the derivation of motion equations is based on a two-module flexible structure. Following 
Lu’s approach, each module of the multi-body structure is divided into several submodules and imaginary beam is 
added in between the center of gravity of two adjacent submodules (see Fig. 7). This idea is, to some extent, 
similar to mass-spring model used in mooring system analysis (Huang, 1992). The difference is that the latter 
ignores the bending and torsion effects and thus springs are added in between mass to model the tension in 
mooring system. However, in the present study, all three effects, i.e. tension, bending and torsion are important 
and should be considered, leading to an imaginary beam (instead of a spring) added in between center of gravity 


of two adjacent submodules to take into account these structural effects. 
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Fig. 7. An illustration of set up for hydroelasticity theory proposed by Lu et al. (2016) 


Then the motion equations for a two-module flexible structure in waves can be established by modifying Eq. 


(7) as follows, 
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where {Fy;} (i=l, 2, ..., 2n) is the force (or moment) acting on the center of gravity of the ith submodule caused 
by the deformation of the imaginary beam. As shown in Fig. 7, the hinge only exists between the n'" and (n+1)" 
submodule. Thus {F,,}+#{0} for k=n andn+1. However, for k being other values, {F,,}={0}. 

Next we will give the expression of {Fs}in Eq.(20). Each side of a beam element i-j has six components of 


displacements and forces (shown in Fig. 8): 
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where ( Fx, Fy; Fz;) are forces in the i” cross section, ( Mx;, My; Mz;) are the corresponding bending moments. It 


should be noted that 7 is the imaginary unit only if it is emphasized. Otherwise, i represents the order number. 
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Fig. 8. Definition of coordinate system and parameters for a beam element 


The relationship between {u};;and {F };jis: 
VE es =|K]., {u}, (22) 
where [K];; is stiffness matrix, the form of which is given in the Appendix B. A detailed explanation of the 


stiffness matrix can be found in the reference (McGuire et al. 2000). The size of the matrix is 12x12. According to 
Eq. (22), we can establish the relationship between force and displacement for both beam k-1 and k (Fig. 7), as 
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follows, 


where the superscript represents the order number of the beam and the subscript represents the order number of 


the submodule’s center. [K E and [K los 


etc. are the partitioned matrix of the stiffness matrix shown in 
the Appendix B. It should be note that fu} =fu}. So we may neglect the superscript of the displacement 
matrix. Thus we have the expression {F Sk \ as follows, 


(F =-{P CP -ia ET Khal [a tol a TL, 25) 


By solving Eq. (20), we obtain the displacement and force of the center of gravity of each submodule. Then 


the displacement {u}; and force {F}; at any position of the structure (suppose that the position is on the K” beam) 


can be calculated as follows, 


P =[K]; ({F}, -[K], e} ) (26) 


{F}, =[K], 4}, HELK] (Fh [Aa 4) 
Eq. (26) is obtained using the similar logic indicated in Eq. (24). The difference is that the variables (i.e. 
displacement, force and stiffness matrix) on the right side of the k” imaginary beam (represented by the order 
number k+1) in Eq. (24) is replaced by the variables on a given position along the k” beam (represented by the 
order number /). The validation of the present approach for a hinged two-module flexible structure (the model 
adopted in Fu et al., 2007) can be referred to Xuet al. (2017). 


3.2. Extension of the approach for flexible structure with arbitrary shape of cross section 

In section 3.1, the flexible structure has a rectangular cross section and the geometric feature of the cross 
section is unchanged along the longitudinal direction of the beam. Therefore, there is an analytical stiffness matrix 
for such a simple beam element (which is given in Appendix B). However, in reality, VLFSs such as floating 
bridge or airport, have very complicated geometric features of cross section. And this geometric feature may 
change along the longitudinal direction of the structure. Thus to extend the hydroelasticity theory shown in 
Section 3.1 to be applicable for multi-module flexible structures with arbitrary shape of cross section is of great 
significance to engineering applications. In the following, we will extend the present approach to be applicable for 
multi-module flexible structures with arbitrary shape of cross section. The geometric feature of the cross section 
may change along the longitudinal direction of the beam (see Fig. 9). The extension of the approach actually 
includes two steps: (1) to represent the physical variables (displacement and forces) of all the nodes of the beam 
element by the physical variables of nodes on the cross section of the beam element; and (2) to represent the 
physical variables of nodes on the cross section by the ones at the center of cross section of the beam element. 


A beam element 


center 


center j 


Fig. 9. A beam element with arbitrary shape of cross section and changed geometric features along the longitudinal direction. 


For a beam element with arbitrary shape of cross section and changed geometric features along the 


longitudinal direction, finite element method (FEM) is applied to the space structure for calculation of the overall 
stiffness of the structure. Suppose that the beam element is discretised into a number of elements and the number 


of nodes is nı. Each node has three components of displacement and force, which is defined as {u*}; = 
(uži Uy; Uzi T and {F*}; = (Fi; Fii Fi) (i = 1, 2, ..., nı). The force is related with the displacement by the 


spatial stiffness matrix [K*] as 


[K Tin Wha 7 Flaa Cn 


In Eq. (27), the nodes are spatially distributed within the beam element. This means that there are nodes both 
on the cross section of the beam (the nodes on the boundary) and in the inner space of the beam. So the first step 
is to represent the physical variables (displacement and forces) of all the nodes of the beam element by the 


physical variables of nodes on the cross section of the beam element, i.e. to establish the stiffness matrix of the 
nodes on the cross section, [Kg]. Eq. (27) can be modified as 
[eh el fe 


el celles i 


where {uz} and {u;} are the displacement of the nodes on the cross section and in the inner space of the beam 
element (shown in Fig. 9), respectively. {F$} and {Ff} are the external force of the nodes on the cross section 
and in the inner space of the beam element, respectively. [K*]pp, [K*]g;, [K*]; and [K*],, are the partitioned 
matrix of the matrix [K*]. It should be noted that {F;}= 0 because these nodes are inner nodes of the beam 
element and no external forces are exerted on inner nodes. From Eq. (28), we can obtain the following 
expressions, 

[Ke Ju} = {Fe} 

Se * * xl * 
[K] =[K lhe [R Ja [x I [x Je 


The second step is to represent the physical variables of nodes on the cross section by the ones at the center 


(29) 


of cross section of the beam element. As indicated in Eq. (21), the equivalent displacement and forces on the 

center of two cross sections at the end of the beam is defined as {u},2,, and {F}i2x1; respectively. If the 
equivalent stiffness is defined as [K], then 

[K]{u}={F}) (30) 

The problem now is to seek the relationship between the equivalent stiffness matrix [K] and the stiffness 

matrix of nodes on the cross section, [Kz]. The displacement of nodes on the cross section can be related to the 


displacement at the center of the cross section as follows, 
{us} [CH (31) 
where [C] is a transform matrix. If the nodes on the cross section is m, then the dimension size for {up}, [C] 
T 
and {u} is 3n7x1, 3n.x6 and 6X1, respectively. The matrix [C] can be expressed as [C K({C,] (C5) cs [c,,]) ; 


where [C;] (the dimension is 3x6) is the partitioned matrix component of [C] corresponding to the 
transformation between the displacement of the i” node and the center. Assuming that the relative position of the 


i™ node to the center of the cross sectionis {rġ;} = (x; yi zi)", then the expression of [C;] is given as follows, 


[G]=|0 1 0 -z 0 x (32) 

i i 3x6 
As the force on the nodes of the cross section and the relative coordinate of nodes to the center of cross 
section are both known, it can be easily to equivalate the force on the nodes, {F$}, to the force on the center of 


the cross section, {F} by the following expression, 
{F}=[CT {Fo} (33) 
Combining Eqs. (29)-(33), we have 
[K]=[c} |x; ][c] 
[k ][e E], 


Eq. (34) shows the relationship between the equivalent stiffness at the center of the cross section and the stiffness 


(34) 


matrix of all nodes of a beam element. Using Eq. 34), the dynamic response of an interconnected multi-module 
floating structure with arbitrary geometric features of cross sections can be solved. 

In order to illustrate the above-mentioned methodology, we choose a cubic steel model with the dimension of 
10m (length) x4m (width) x3m (height). The other parameters are given as follows: the Young’s modulus 
E=2.1x10'' pa, the Poisson’s ratio 0.3 and the density 7850kg/m*. Using the analytical expression given in 
Appendix B, the 12x12 stiffness matrix [K ]qna is given in Table 1. The stiffness matrix is also calculated using 
the FEM method, as indicated in Eq. (34). The steel model is discretised with the number of nodes 20, 8 and 6 in 
the direction of length, width and height, respectively. The ANSYS is used to extract the stiffness matrix [Kp] 
and the equivalent stiffness on the center of the cross section [K ]pmy is calculated using Eq.(34). The results are 
given in Table 2. The difference between [K ]ana and [K ]Fem is given in Table 3. It can be seen that the relative 
error of all components of the stiffness matrix calculated by two approaches is within 10%, which validates the 
correctness of Eq. (34). 

Next we consider the model adopted in experiment conducted by Yago and Endo (1996) to further validate the 
above-mentioned methodology. The model is a continuous pontoon-type VLFS with the length, width, height, 
draught and mass being 300, 60, 2, 0.5m and 9.225x10° kg, respectively. The bending rigidity is EI= 4.77x10" 
Nm’. The hydroelastic response of the model is calculated for heading waves of different wave length. The 
structural stiffness is calculated using both analytical expression in Appendix B (denoted as ‘Approach 1”) and FEM 
method with Eq. (34) (denoted as ‘Approach 2”). The experimental results obtained by Yago and Endo (1996) is 
denoted as ‘Experiment’. As shown in Fig. 10, the dynamic response of the continuous VLFS in waves is almost 
identical for two approaches, both of which agrees well with experimental results. This further validates that the 
correctness of extension of the hydroelasticity theory proposed by Lu et al. (2016). 


Table 1 the stiffness matrix [K ]ana calculated by analytical expression 
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Table 2 the stiffness matrix [K ]rem calculated by FEM method (Eq.(34)) 
1 2 4 5 6 7 9 Ped pes 
2. 52E+11 0 0 0 0 
0 0 0 0 1. 35E+11 
0 0 —1. 8E+10 —8. 9E+10 0 
0 0 0 1. 57E+11 0 0 0 0 0 0 
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=2. 5E+11 0 0 0 0 2. 52E+11 0 0 0 
0 -2. TE+10 0 0 —1. 3E+11 0 0 0 =} _ SEEI 
0 0 —1. 8E+10 0 8. 87E+10 0 0 1. 77E+10 8. 87E+10 0 
0 0 -1. 6E+11 0 0 0 0 0 0 
0 0 -8. 9E+10 0 2. S4E+11 0 0 8. 87E+10 6. 32E+11 0 
0 1. 35E+11 0 0 3. 37E+11 0 0 0 1. 01E+12 
Table 3 The difference between [K]ana and [K ]rem (%) 
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0 0 5. 66 0 10.3 0 0 0 5. 66 3.9 0 
0 8. 47 0 0 0 18. 48 0 8. 47 0 0 5.5 
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Fig. 10. Proof of correctness of extension of hydroelasticity theory proposed by Lu et al. (2016) 


4. Time-domain analysis of hinged multi-body structure with moving point loads on upper 
surface 


4.1. Motion equationof hinged two-module flexible structure in time domain 
Frequency domain method mentioned above is invalid if nonlinear external loads or transient loads need to 
be considered. In such conditions, the hybrid frequency/time domain method ((i.e. the time domain method based 
on impulse response function (IRF)), which was first introduced by Cummins (1962), is adopted to establish the 
motion equations. Here we consider a hinged two-module flexible structure in waves. The motion equations of the 
structure (which is divided into 2n submodules) is established using the SAM approaches as follows, 
[mlao] [ato] = [nae] tix) 
[an (<0) | [M,]+[ a. (20) | oy [asan (x) | {iiO} 


[azna (<0) | ca (<0) | Si | Monon |+[ aonan (<0) | {ü} (35) 
[Bi(t-7)] [Bolt-2)] e [Bato] nO fa} Ma 
j [2u] [Balt] a [Balt] (a0) dr +Kan] t Mea) 


12nx12n 
0 : 


[Ba (t-1)] [Baat] e [Banan t-a] | tn} funt {Foon} 
where [a; j (co) is the added mass (or moment) matrix of the i" module caused by the motion of the j” module at 
frequency of infinity. [B; KO) is the kernel function of the i'" module caused by the motion of the jf” module, 
which is related to the radiation damping matrix, [b; j(w)]. The expressions of and la; j (co) | and [B; j O) are 


given as follows, 


[a, (o)]= [a (2) | [B OJo, nat 


(36) 


2 po 
|B, (t) | = =! |b, (@) |cos(@t)de 


where œ, is an arbitrarily chosen frequency. 


The total stiffness matrix [Krotaı] is composed of three parts, i.e. the hydrostatic stiffness of the structure [C], the 
structural stiffness matrix (considering the deformation) [K] and the approximate stiffness matrix (using SAM 
approach to approximate the hinge constraint of the structure) [Kcons]. The expression of [Krotaı] is given as 


follows, 


[Kroi larza =[C] T [x] w [Koons ] 
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(37) 


[Oo] 


The model adopted by Fu et al. (2007), which is a hinged two-module flexible structure, is used here to 
validate the time domain approach. The hinged two module structure is based on a continuous pontoon-type VLFS 
with the main parameters given in Section 3.2. Fu et al. (2007) separated the continuous model into two modules 
and added a hinge connection in between two modules. The dynamic response of the hinged two-module flexible 
structure in waves of different wave length using both frequency and time domain method is given in Fig.11. It 
can be seen that the agreement between frequency domain and time domain is quite good, which validates the 


correctness of time domain approach shown in Eqs. (35)-(37). 
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Fig. 11. Frequency (FD) and time (TD) domain analysis of dynamic response of a hinged two-module flexible structures in waves. 


4.2. Consideration of a moving point load 

For a floating bridge or airport, there are vehicles moving on the upper surface of the structure. Usually the 
size of vehicles are small compared with the floating structure, which means that the vehicle can be regarded as a 
point load. Thus the problem of a vehicle moving on the floating bridge (or airport) can be simplified as a point 
load moving on the upper surface of a hinged two-module flexible structure. In this section, the time domain 
method is adopted to investigate the dynamic response of a hinged two-module flexible structure with moving 
point loads on the upper surface in waves (see Fig. 12). 

A point load 


U, 
— > 


—> 


Incident waves A hinged two-module flexible structure 


Fig. 12. A schematic of a hinged two-module flexible structure with moving point mass on the upper surface in waves 


The linear assumption for the wave structure interaction still remains here. Besides, we assume that the 
existence of the moving point mass does not cause large structural deformation. Thus the floating system shown in 
Fig. 12. is still a linear system. The point mass is assumed to move in a uniform horizontal velocity U,. The force 


exerted by the point mass on the structure is denoted as {hy Ohr As the motion equation is established on the 


center of each submodule of the structure. Thus the force { f O} should be equivalated to the ones on the center 


of submodules, which is denoted as { ho Ohana =C foi O { h2 (t)} ... { em), J: 


In order to obtain the relationship between { f (t)} and { to O}, we first consider a (static) point load acting 
on a hinged two-module flexible structures in calm water. The motion equation (can be obtained from Eq. (35) 


after some manipulation) is given as follows, 


Ver ea} 38) 


A ba) 


It should be noted that in calm water, each module of structure is divided into n; submodules. n, is a large 


[ Total les, x12n, 


number (n;=150 adopted in this section) and is much larger than the number of submodule used in hydroelastic 
reponse,  (n=4 in the present analysis). The reason is that by setting a series of submodules (or elements) of very 


small size, the point load can be regarded as acting on the center of one particular submodule (for example, the 
if” submodule shown in Fig. 13). Then the equivalent force satisfies, es j}=0 (j#i,) and os gah (t)} 
(=i). 


A point load 
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A hinged two-module flexible structure 


Fig. 13. A schematic of a static point load acting on a hinged two-module flexible structure in calm water 


By solving Eq.(38), we obtain the distribution of displacement along the structure in calm water caused by 
the static point load at a given position on the structure. Then the displacement at the center of submodules (the 


number is 2n), {U}42nx1 for hydroelastic calculation (see Eq.(35)) is also obtained. Thus we can obtain the 


equivalent force { i} as follows, 


12nx1 


(39) 


a [ Kea are 


T (un) 


In calm water, we can obtain the equivalent force { fo} i for static point loads acting on any position of 


12nx 


the structure and establish a database to store these values. When we establish the database, the point load is 
chosen as a unit load. Then the strategy of dealing with dynamic response of a hinged two-module flexible 
structure with moving point loads is summarized in Fig.14. The point load considered is a vertical load, which is 
the sum of the inertial force and gravity force. ü, represents the vertical acceleration of the structure at the 
position where the point load is exerted. Here we assume that the floating system is a linear system, which means 


that the superposition principle is applicable. Thus for every time step, we can add the equivalent force caused by 


the static point loads in calm water to the time domain equation of the structure in waves. 


1. Establish a data base to store the equivalent force { founi) for a unit point 


12nx1 
load acting on any position of the hinged two-module flexible structure in calm 


water, using Eq. (38) and Eq.(39). 


2. At time t, the point load is Q ,=m(ü; — g). The position of the point load is x,. Then extract the 
force { founi (from the database. The equivalent force { fo (tn)} =Q,,X { Founit} oma’ 


12nx1 


, 10 the right hand side of Eq.(35) and solve the time 


12nx 


3. Add the equivalent force {f, (ta) } onc 


domain equation for next time step fp. 


4. Repeat step 2 and 3 until the time instant t > T (T is the duration of simulation). 


Fig. 14. The strategy of dealing with dynamic response of a hinged two-module flexible structure with moving point loads 


A continuous flexible VLFS model with the same parameters adopted in Yago and Endo (1996) (see Section 
3.2) is adopted to validate the above-mentioned numerical procedure. The only exception is that the length of the 
model is changed to be 1000m, representative of an infinitely long model. We consider a static point load of IMN 
on the middle of the VLFS model in calm water. This problem is actually equivalent to a static load exerting on a 
beam of infinite length and on elastic foundations (The hydrostatic restoring force due to change of draft acts as 
an elastic foundation). And the latter is a classic problem in structural mechanics and analytical solution is 
available (Timoshenko, 1940). As the hinge connection is not considered, the third term in the total stiffness 
matrix caused by the hinge (see Eq. (37)) should be cancelled. This problem is not directly relevant to the scope of 
this section, but it is treated as a preliminary test case to validate the numerical procedure for treatment of moving 
loads on hinged multi-module structures. In our calculation, the continuous model is uniformly divided into 200 
modules. The numerical results are compared with analytical ones, which (in Fig.15) is labelled as “Present” and 
“Theory”, respectively. It can be seen from Fig.15 that the agreement is quite good, which partially validates the 


numerical procedure. 
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Fig. 15. Results of a static point load acting on an infinitely long VLFS model: (a) vertical displacement; and (b) bending moment. 


Next we consider the moving load experiment conducted by Endo and Yago (1999). The experimental modal 


structure, VL-10, was so designed as similar to the continuous flexible VLFS model (MF-300) adopted in Yago 


and Endo (1996) on the basis of similitude (with a scaling factor of 30.77). In the following description of 
experiment, all physical properties are listed in prototype scale (i.e. the scale of MF-300). In the experiment, a 
movable weight of 201016kg was translated on the runway from x/L=-0.31~0.4 (see Fig. 16) with almost constant 
velocity, 3.39m/s. It should be noted that the origin is set in the middle of the structure in Endo and Yago (1999), 
which is different from the definition of origin at left side of the structure in previous sections in the present study. 
This case is also calculated by the present numerical method. It should be noted that the third term in the total 
stiffness matrix caused by the hinge (see Eq. (37)) should be cancelled as this is a continuous structure. Besides, 
as no incident waves are considered, the wave excitation force in the right hand side of Eq. (35) should be ignored. 


The time history of vertical displacement of the structure obtained from the present numerical method, “Present”, 


and from the experiment, “Experiment”, are comparatively shown in Fig. 17. It can be seen that the calculated 


results are found to show a reasonable agreement with the experiment, which further validates the proposed 


numerical methodologies for calculating the dynamic response of a flexible structure undergoing moving point 
loads. 
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Fig. 17. Time history of vertical displacement of a continuous flexible structure excited by a moving load in calm water 


After validation of the present numerical model, we consider a hinged two-module flexible structure 
undergoing a moving point load in waves. The same VLFS model adopted in Yago and Endo (1996) is used here 
(see Fig.16). The continuous VLFS is divided into two identical modules and a hinge connection is added in 
middle of the VLFS. Here, for consistency, we still set the origin of the coordinate system at the left side of the 
structure. A movable mass of 800000kg was translated on the left side of the structure with constant velocity, 5m/s 
from 40s. The amplitude of incident wave is 1m. The wave length is 420m (the period is 19.55s for water depth of 
58.5m considered here). The vertical displacement of the structure excited by the combination of wave and 
moving load is given in Fig.18. It can be seen from Fig. 18 (a) that before the point load is exerted on the structure, 
the deflection (or the vertical displacement) of the structure is periodic as the incident wave is regular and periodic. 
The contour of response of the structure is characterized by periodically varied oblique band of different colour. It 
can be seen that the maximum structural deflection occurs at the middle (where the hinge is applied) and the ends 
of structure (where the deepest red colour of contour exists in Fig.18 (a)). From 40s, a point load is exerted on the 
structure and moves forwards with a uniform velocity. The dynamic response of the structure is affected by the 
moving load, with positive (upward) deflection decreased and negative (downward) deflection increased. The time 
series of structural deflection at different locations of the structure together with input waves are shown in Fig.18 
(b)-(g), which shows clearly the effects of the moving point load on the dynamic response of the structure. It 
should be noted that in this section, we do not intend to provide an extensive analysis of the dynamic response of a 
flexible structure with moving point loads for different parameters of wave, structure and point loads. Our purpose 
here is just to choose a case to validate the effectiveness of the proposed strategy for dealing with such a complex 
problem. 

The bending moment of the structure excited by the combination of wave and moving load is given in Fig. 19. 
The contour of the bending moment (Fig.19 (a)) highlighted the significant effects of moving point loads. The 
maximum bending moment is almost doubled due to the influence of moving point loads (compared with the 
bending moment of structure induced only by waves). Fig.19 (b)-(d) show the time series of bending moment for 
different locations along the longitudinal direction of the structure. As both ends of the structure is free, the 
bending moment at x/L=0, 0.5 and 1 is 0, which is not shown in Fig. 19. It can be seen that before the point load is 
put on the structure, the time series of bending moment is smoothly sinusoidal. After the point load is exerted, an 
increase of bending moment occurs (at x/L=0.25, the bending moment is almost doubled). And high-frequency 
oscillation also appears as a disturbance to the wave-frequency time series of bending moment. To the authors’ 
knowledge, no experimental results for such cases are available in open literature. Therefore there is no 


comparison between the present numerical results and other experimental data. 


5. Conclusions 


In this paper, we developed a consistent approach which is capable of calculating the dynamic response of 
rigid/flexible multi-body floating systems. This method is based on multi-rigid-body hydrodynamics and 
Euler-beam bending theory and can be used to deal with rigid or flexible multi-body floating systems in both 
frequency and time domain. This method is consistent in a sense that the dynamic response of rigid/flexible 


structure can be calculated by just turning off/on a structural stiffness matrix in in multi-rigid-body motion 


equations. Unlike traditional modal expansion approach, the present method does not require the calculation of 
flexible motion mode and it is more computationally efficient. The present method can be applied more easily in 
industry practice with the aid of commercial code of hydrodynamics such as HY DROSTAR. 

Within the framework of this consistent approach, a stiffness approximation method is proposed to calculate 
the natural frequency and corresponding oscillation mode of a hinged multi-module structure. Instead of directly 
describing the constraint condition of a hinge connection, a stiffness of very large value is adopted to approximate 
the constraint condition. The determination of natural frequency of a multi-body floating system is significant for 
the purpose of structural safety of traditional multi-body floating system or power efficiency of wave energy 
converters. 

The hydroelasticity theory proposed by Lu et al. (2016), which is used to deal with continuous flexible 
structures with cross section of simple shapes, is extended by using the ideal of finite element method, to be 
applicable for multi-module flexible structures with cross section of arbitrary shapes and varied geometric features 
in longitudinal direction. Validation of the extension has been performed. The extended hydroelasticity theory 
makes it possible to calculate the dynamic response of more complex flexible structures such as floating bridges 
and containership. 

Finally, an effective numerical strategy is established to analyse in time domain the dynamic response of a 
hinged multi-module flexible structure with moving point loads on the upper surface in waves, representative of 
vehicles moving on a floating bridge or airport. Validation of the numerical strategy has been made for a static or 
moving point load exerted on a continuous flexible structure in calm water, which is a less complicated problem. 
The structural deflection and bending moment of a hinged two-module flexible structure undergoing both wave 
and moving point loads calculated by the present numerical approach are given, which highlighted the effects of 
moving point loads. In future work, relevant model testing will be carried out to further validate the present 


numerical approach. 
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Fig. 18. The vertical displacement of the structure excited by the combination of wave and moving load. (a) the contour of vertical 


blue solid line represents structural deflection and the red dash line indicates the incident wave. 


displacement with respect to time; 


— - Trajectory of moving load 
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Fig. 19. The bending moment of the structure excited by the combination of wave and moving load. (a) the contour of bending 


moment with respect to time; (b)-(c) are time series of bending moment at x/L= 0.1, 0.25, 0.4, respectively. 
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AppendixA Proof of Eq. (11). 


For a multi-body floating system without constraints, the motion equation is established as follows (similar 
as Eq. (8)), 

[K](u} ={F} (Al) 
where [K], {u} and {F} are the overall stiffness matrix, displacement vector and external force, respectively. 
For a floating system of m modules, the dimension size for [K], {u} and {F} is 6m x 6m, 6mx1 and 6mx1, 
respectively. 

The total energy II, which is the sum of the elastic strain energy (stored in the deformed body) and work 
done by external force, is given as follows, 
T= fu)" [uh fo)" {F} (A2) 


The motion of the floating system is constrained by a constraint matrix [L] as follows, 


[L] tu} =0 (A3) 
According to the minimum total potential energy principle, the solution of the displacement {u}, is 


equivalent to the mathematical problem: to find the extreme value of II with the constraint indicated in Eq. (A3). 


Then the method of Lagrange multiplier can be applied. We introduce a Lagrange multiplier {Ẹ} and then define 


the Lagrange function IT as follows, 


I= stay" [k] te) LF} fo} [Ll (A4) 
Then the variation of Eq. (A4) gives 
oT = {su} ([K]u} -{F} -[2]' {}) {00} L (AS) 


Therefore, we have 


(A6) 


T 

[K]iu}-[Z] te} =1F} 
[E}tu} =0 

By comparing Eq. (A6) with Eq. (12), we find that {ọ} is actually the force on the connection (or constraint 

element) {fz}. And by comparing Eq. (A6) with Eq. (7), we have 


{F} =I] {oE {fc} (A7) 
The proof of Eq. (11) is finished. 


AppendixB Stiffness matrix for a uniform beam element 


The stiffness matrix [K].; in Eq. (21): 
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The size of the symmetric matrix is 12x12. The matrix is affected by the size of each submodule, Young 
modulus and Poisson ratio. / is the length of beam element (the dimension along x axis). b is the dimension along 


y axis. h is the dimension along z axis, «is a constant related to b and h. A is cross-sectional area. E is Young 


modulus. Jis the correction factor considering shear deformation for the beam element bending in xoy plane. £, is 
the correction factor considering shear deformation for the beam element bending in xoz plane. /, is the moment of 


inertia around the y axis. /, is the moment of inertia around the z axis. G is shear modulus. 
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